## Base RF model
set.seed(42)
library(ggplot2) # for autoplot() generic
library(dplyr)
library(ggpubr)
library(ggcorrplot)
library(sf)
library(stars)
library(tmap)
library(car)
dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 31 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS: GCS_unknown
my_vars <- c("mb_mwea", "area_km2", "z_med", "hi", "z_aspct", "z_slope", "tau",
"t2m_18", "t2m_d", "tp_18", "tp_d", "tpseas_18", "tpseas_d")
my_labels <- c("Mass balance", "Area (km2)", "Median elevation", "Hypsometry",
"Aspect", "Slope", "tau",
"2m temperature (2018)",
"2m temperature change (2018 - 2008)",
"Annual precipitation (2018)",
"Annual precipitation change (2018 - 2008)",
"Precipitation seasonality (2018)",
"Precipitation seasonality change (2018 - 2008)"
)
n_vars <- length(my_vars)
dat <- dat_sf %>%
st_drop_geometry() %>%
dplyr::select(any_of(my_vars))
dat_cor <- cor(dat)
ggcorrplot(dat_cor)
Variance inflation factors (\(>5\) is a conservative threshold for multicollinearity).
fit <- lm(mb_mwea ~ ., dat)
vif(fit)
## area_km2 z_med hi z_aspct z_slope tau t2m_18 t2m_d
## 1.344101 1.760093 1.098699 1.009871 1.348702 1.928942 2.727786 1.387143
## tp_18 tp_d tpseas_18 tpseas_d
## 2.330231 1.721179 2.423583 1.483386
for(i in 1:n_vars) {
p1 <- gghistogram(dat, my_vars[i], main = my_labels[i])
print(p1)
}
Possible variables for log-transform:
area_km2z_slopetautp_18dat$larea_km2 <- log(dat$area_km2)
p1 <- gghistogram(dat, "larea_km2")
print(p1)
dat$lz_slope <- log(dat$z_slope)
p1 <- gghistogram(dat, "lz_slope")
print(p1)
dat$ltau <- log(dat$tau)
p1 <- gghistogram(dat, "ltau")
print(p1)
dat$ltp_18 <- log(dat$tp_18)
p1 <- gghistogram(dat, "ltp_18")
print(p1)
Load basemap
basemap <- read_stars("~/Dropbox/Data/summer/natural_earth/HYP_50M_SR_W/HYP_50M_SR_W.tif")
basemap2 <- st_crop(basemap, st_bbox(dat_sf))
basemap2 <- st_rgb(basemap2)
for(i in 1:n_vars) {
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(dat_sf) +
tm_symbols(col = my_vars[i],
size = 0.25,
alpha = 0.75,
style = "quantile") +
tm_layout(main.title = my_labels[i],
legend.position = c("left", "bottom"),
legend.bg.color = "white",
legend.bg.alpha = 0.5)
print(m1)
}